import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
np.random.seed(908) # a not-so-random starting seed
Let's start with a uniform distribution on the unit square [0,1]×[0,1] . Create a 2D array samples of shape (2, nsamples):
nsamples = 10**2
samples = np.random.rand(2, nsamples)
Scale the sample points "samples", so that we have a uniform distribution inside a square $[-1,1]\times [-1,1]$. Calculate the distance from each sample point to the origin $(0,0)$
xy = samples * 2 - 1.0 # scale sample points
r = np.sqrt(xy[0, :]**2 + xy[1, :]**2) # calculate radius
## talk about here how we can reduce the calculation by removing the sqrt
plt.plot(xy[0,:], xy[1,:], 'k.')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y');
If the distance is less than 1 (radius of the circle), then the sample point is inside the circle.
incircle = (r <= 1.0)
n = np.arange(1, nsamples+1)
countincircle = 4 * incircle.cumsum() / n
The approximation for $\pi$ when we use n
sample points is:
pi_approx = countincircle[-1]
print(pi_approx)
plt.plot(xy[0,np.where(incircle)[0]], xy[1,np.where(incircle)[0]], 'b.')
plt.plot(xy[0,np.where(incircle==False)[0]], xy[1,np.where(incircle==False)[0]], 'r.')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y');
n
(number of sample points)¶plt.plot(n, countincircle, '.')
plt.xlabel('n')
plt.ylabel('count in circle');
error = np.abs(countincircle - np.pi)
plt.loglog(n, error, '.')
plt.xlabel('n')
plt.ylabel('error')